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Abstract 

The reconstruction of large phylogenetic trees from data that violates 
clocklike evolution (or as a supertree constructed from any m input trees) 
raises a difficult question for biologists - how can one assign relative dates 
to the vertices of the tree? In this paper we investigate this problem, as- 
suming a uniform distribution on the order of the inner vertices of the tree 
(which includes, but is more general than, the popular Yule distribution 
on trees). We derive fast algorithms for computing the probability that 
(i) any given vertex in the tree was the j-ih speciation event (for each 
j), and (ii) any one given vertex is earlier in the tree than a second given 
vertex. We show how the first algorithm can be used to calculate the 
expected length of any given interior edge in any given tree that has been 
generated under either a constant-rate speciation model, or the coalescent 
model. 

Keywords: Phylogenetics, neutral model, dating speciation events, edge 
lengths. 

1 Introduction 

A fundamental task in evolutionary biology is constructing evolutionary trees 
from a variety of data. These constructed trees show the ancesteral relationship 
between the species. 
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Not only the relationship between species is of interest, but also the time 
between speciation events. When constructing an evolutionary tree from a set 
of molecular data which satisfies the molecular clock, the edge lengths can be 
interpreted as a time scale. In many cases, no time scale is obtained when 
constructing a tree though: 

• Often, molecular data does not satisfy the molecular clock and so the edge 
lengths do not represent a time scale. 

• Trees can be constructed from morphological data or non-standard molec- 
ular data like gene order. This does not provide any edge lengths. 

• Having several different trees, one can combine them and construct a 
'supertree'. Even though there may have been time scales on the original 
trees, most supertree methods return a tree without a time scale. 

For those trees, we still want to find edge lengths representing the time 
between speciation events. In this paper, we will estimate the edge lengths from 
the shape o f the tree. The method works for trees which evolve d under the 
Yule model \Yu\(l ll924L fEdwardsL Il97nl iHardind. Il97ll IPatrel IT99ll . Under the 



Yule model, in each point of time, each species is equally likely to split. Minor 
changes to the method for the Yule model give us an edge le ngth estimation for 
trees under the popular coalescent setting |Nordbor'3 l200l| . 

An example for a tre e with unknown e dge lengths is the primate supertree 
Tp recently published in |Vos and Mooers| . Figure ^ shows a part of Tp. The 
primate tree is a supertree on 218 species and was constructed with the MRP 
method (Ma trix Representation using Parsimony analysis, see Bauml Il992l 



iRaeanL ri992l |). Since for most of the interior vertices, no m olecular estimates 
were available, the edge lengths for the tree were estimated. In |Vos and Mooersl| . 
10^ rank functions on Tp were drawn uniformly at random. For each of those 
rank functions, the expected time intervals, i.e. the edge lengths, between ver- 
tices were considered (the expected wa iting time after th e (n — l)th event until 
the nth event is 1/n). The authors of |Vos and Mooersj concluded their paper 
by asking for an analytical approach to the estimation of the edge length, which 
we will provide below. 

In order to estimate the edge lengths, we developed the algorithms RankProb 
and Compare. Those algorithms answer questions like: 

Was speciation event with label 76 in the primate tree (see Fig. ^ more 
likely to be an early event in the tree or a late event? What is the probability 
that 76 was the 6th speciation event? Was it more likely that speciation event 
76 happened before speciation event 162 or 162 before 76? 

The algorithms work for trees where every labeled history is equiproba- 
ble. This class of model, which includes the Yule mod el and the coa, l escent 
model, has been popular i n mac roevolutionary studies Nee and MavL Il997i 



[zhaxybavcva and Gogartc^, l2004j . Note that the algorithms here are the same 
for the Yule model and the coalescent model, whereas the edge length estimation 
has minor differences for the two models. 
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Figure 1: Part of the prima te supertree. Figure 4-13 are some subtrees, for 
details see IVos and Mooersl. 



The algorithms RankProb, Compare and an algorithm for obtaining the 
expected rank a nd variance for a vertex were implemented in Python, see 



2 Probability distribution of the rank of a ver- 
tex 



Let T be a rooted phylogenetic tree [Semple and SteelL l2003| with \V\ — 



n 



leaves. The set of interior vertices of T shall be V . For a binary tree, we have 
|F| = n — 1. Let the function r be a bijection from the set of interior vertices 
y of T into {1, 2, . . . , \V\] with r{vi) < r{v2) if vi is an ancestor of V2- The 
function r is called a rank function for T. A vertex v with r(v) = z is said to 
have rank i. Note that r induces a linear order on the set V. Further, define 
r{T) :— {r : r is a rank function on T}. We are interested in the distribution 
of the possible ranks for a certain vertex, i.e. we want to know the probability 
of r{v) = i for a given v d V. If every rank function on a given tree is equally 
likely, we have 

\nT)\ 
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Figure 2: Labeling the tree for the algorithm RankProb 



which will be calculated for rooted binary trees in polynomial time by algorithm 
RankProb. In the algorithm, we will use the formula Semple and Steel . .2003j 

where Uv is the number of leaves below v. Note that Equation[21holds for binary 
and nonbinary trees. 

Examples of stochastic models on phylogenetic trees where each rank function 
is equally likely include: 

• The Yule model has the probability distribution P[r|T] = 
which is the uniform distribution lEdwardfl FrnTfi iRrowuL [^99^ . 

• The coalescent model has the same probability distribution on rooted bi- 
nary ranked tr ees as the Yule model. So P[r|T] is the miiform distribution 
|Aldouil200l| . 



For s ome sets of trees (e.g. those drawn from the uniform model Pinelisl 



|2003| . also known as PDA model), no rank function is induced. If one 
assumes that all rank functions are equally likely on these trees, one can 
apply Equation n to such trees as well. 

2.1 A polynomial-time algorithm 

The following algorithm calculates the probability distribution of the rank of a 
vertex w in a rooted binary phylogenetic tree T. The idea of the algorithm is the 
following (c/. Figure [2Jl. Label the vertices on the path from v to the root p by 
V = xi, . . . ,Xn = p- Let %n be the subtree of T containing the vertex Xm and 
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all its descendants. Let ar„,t)(j) be the number of rank functions on the tree %n 
where v has rank i. The values aT„,t)(*), i = 1, . . . , \V\ are calculated iteratively 
for m = 1, . . . , n. The probability P[r(w) = i] equals — "j^"'^'*-'-* — . The a-values 

in the fraction have a lot of factors in common which cancel out. In the fol- 
lowing algorithm, we calculate a-values without the unnecessary terms instead, 
ar^Ai)- We have ar„,.(«) = aT^,,{i)\r{Ti)\\r{ri)\\r{ri)\ . . . \r{%'^_^)\. 

Algorithm: RankProb(T, ti) 

Input: A rooted binary phylogenetic tree T and an interior vertex v. 
Output: The probabilities P[r(t;) = i] for i = 1, . . . , 
1: Denote the vertices of the path from v to root p with 

{V = Xi,X2, ...,Xn=p). 

2: Denote the subtree of T, consisting of root Xm and all its descendants, by 
7m for m = 1, . . . , n. 

Initialize aT„,t)(*) for i — 1, . . . , |V7-|,r7i — 1, . . . ,n 



ari,«(l) := 1 

for TO = 2, . . . , 77 do 

T^_i %n \ (7;„-i U Xm) {cf. FigureEl 
for i — m, . . . ,\ V-Trr, I do 
M :=min{|l/r'r M " 2} 



i-l I 

M 



end for 
end for 

for i = 1, . . . ,\Vr \ do 
end for 

RETURN P[r(w) =i],i^ 1,...,\V\. 



\Vr'\-J /V ^ 



Proving the correctness and runtime of RankProb makes use of the follow- 
ing two observations. 

Remark 1. Let Ai be a set containing rii elements with a linear order, i G {1, 2}. 
There are ("^,^^"^) possible linear orders on Ai U A2 which preserve the linear 
order on Ai and A2 . This follows from the observation that the number of such 
linear orders on Ai U A2 is equivalent to the number of ways of choosing ni 
elements from ni + 112 elements, which is ("^^j"^)- 

Remark 2. The values (^) for all n,k < N {n,k,N ^ N) can be calculated in 
0{N'^) using Pascal's Triangle. Thus, after 0{N'^) calculations, any value (^) 
with n,k < N can be obtained in constant time. 

Theorem 3. RankProb returns the quantities 

¥[r{v) = i] 
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Figure 3: Labeling the tree for the recursion in RankProb 



for each given v & V and all i E 1, • • ■ , The runtime is 0{\V\'^). 

Proof. Let ar„.„(i) = ar,„,viiMTi)\\r{T{)\\r{T^)\ . . . \r{T;^_^)\. We first show 
that ar„,t)(*) = \{r : r{v) ^ i, r E r{%n)}\ for m = 1, . . . , n, i = 1, . . . , \Vt\- 
That implies 

pr / N -1 |{r : r{v) =i,r e r{T)}\ ^ aT^vji) ^ aT^vji) 
^ \r{T)\ E^c^tA^) E,«r,„W 

which proves the theorem. 

The proof is by induction over m. 
For m = 1, ari,.(l) = |r(Ti)|5ri,. W = HT,)\ = \{r : r{v) - 1, r e r{T)}\. 
Vertex v is the root of 71, so aTi,v{i) = for alH > 1. 

Let m = k and aT„,,i>(i) = \{r '■ r{v) = i, r € r{%n)}\ holds for all m < k. 
OiTk,v{i) — clearly holds for all i > IVtiJ since rg; : v — > {1, . . . , |^Tfc|}- So it 
remains to verify that the term (*) returns the right values for axk.vi'i)- Assume 
that the vertex v is in the {i ~ j ~ l)-th position in Tk-i (with i — j — 1 > 0) 
for some rank function rT^^i and v shall be in the i-th position in Tfe. 

Now combine the linear order in the tree Tk-i induced by rT^^i with a linear 
order in Tl_-^ induced by rq-^ ^ to get a linear order on Tk- The first j vertices 
of Tl_^ must be inserted between vertices of T^-i with lower rank than v so 
that V ends up to be in the i-th position of the tree 7^. Count the number of 
possible way to do this as follows. The tree 7^'_i has |r(7^'_-^)| possible rank 
functions. Combining a rank function r7-^_j^ with a rank function rg-^' ^ to get 
a rank function r-r^ with rq-^{v) = i means inserting the first j vertices of "T^Li 
anywhere between the first {i j ~ 2) vertices of 7fe-i. There are 

possibilities according to Remark ^ For combining the |VTfc_i| ~ (* ^ j ^ 1) 
vertices of rank bigger than v in 7fe„i with the remaining | Vj^ ^\— j vertices in 
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7^'_;^, there are 

|Vt.-x 1 - - J - 1) + \Vri_, I - J\ /|Vr._, I + \Vr^_, 
\Vr-J-J J \ 

possibilities. This foUows again from Remark ^ The number of rank functions 
fTk-i with TTk-iiv) = i — j — I is aTk-i,v(i — j — by the induction assumption. 
Multiplying all those possibilities gives 



l^r,_ J + |^r,L J - - 1)\ A - 2' 
\Vt'\-J J\ J 



where aT._„.(*) = 5r,_i,.(*)IK^i)IIK^i')IIK^2 )l ■ • ■ IK^^z)!- The value \{r : 
r{v) = i, r € r{T)}\ is then the sum over all possible j which establishes the 
correctness of the algorithm. 

All that remains is to verify the runtime. Note that the combinatorial factors 
(^) for all n,k < \V\ can be calculated in advance in quadratic time, see Remark 
13 In the algorithm, those factors can then be obtained in constant time. 

The most time consuming part of the algorithm is line 13. Adding up all 
calculations needed for obtaining a'j- ^(i), m = 1, . . . , n, i = 1, . . . , \Vt^ \ comes 
to: 



m— 2 m— 2 m— 2 

The last inequality holds since the vertices of the , m=l....,n— 1, are 
distinct. Therefore, the runtime is quadratic. □ 

Remark 4. With P[r(f;) = i] from Theorem 13 the expected value fJ,r(v) and 
the variance cr^^^^ for r{v) can be calculated by 

Remark 5. The algorithm RankProb can be generalized to non-binary trees 
GernhardlEif)^ . The runtime is again quadratic. 



3 Application of RankProb - Estimating edge 
lengths 

3.1 The Yule model 

A very common stochastic model for rooted binary phylogen etic trees with edge 
lengths is the continuous- time Yule model [Edwardi. ll97Clj . As in the discrete 
Yule model, at every point in time, each species is equally likely to split and 
give birth to two new species. The expected waiting time for the next speciation 
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Figure 4: Labeling the tree for estimating the edge lengths. 



event in a tree with n leaves is \/n. That is, each species at any given time has 
a constant speciation rate (normalized so that 1 is the expected time until it 
next speciates). 

Assume that the primat e tree Tp evolved under the continuous-time Yule 
model. In GernhardL |2006| . the tree shape of 7^ (i.e. the tree without edge 



lengths) under the discrete Yule model is tested against the uniform model and 
accepts the Yule model. 

Here, we describe how to estimate the edge lengths for a tree which is as- 
sumed to have evolved under the continuous-time Yule model. 

Let {u,v) be an interior edge in T with u the immediate ancestor of v. Let 
X be the random variable 'length of the edge (ii,w)' given that T is generated 
according to the continuous-time Yule model. 

The expected length E[X] of the edge (u, v) is given by 

E[X] ^^¥\X\r{u) = i,r{v) = j]F[r{u) = i,r{v) = j]. 

Since, under the continuous-time Yule model, the expected waiting time for the 
next speciation event is 1/n it follows that: 

fe=i 

It remains to calculate the probability P[r(u) = i,r{v) ~ j]. This is equivalent 
to counting all the possible rank functions where r(u) = i and r[v) = j. The 
subtree % consists of v and all its descendants. The tree 7^ equals the tree T 
where all the descendants of v are deleted, i.e. w is a leaf in 7^, see Fig. 0] 

Note that V[r{u) = i,r{v) = j] = if \Vt^\ < j — 1. Therefore, assume 
I^Tu I > i — 1 in the following. 

The number of rank functions on is |r(7^)|. The probability P[r(u) = i] 
can be calculated with RankProb(7^, u). So the number of rank functions in 
Tu with V[r{u) = i] is P[r(u) = i] ■ |r(r„)|. 

The number of rank functions on Ty is |r(7^)|. Let any linear order on the 
trees 7^ and 7^ be given. Combining those two linear orders into an order, r, 
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on T with r{v) = j means that the vertices with rank 1, 2, . . . , j — 1 in 7^ keep 
their rank. Vertex v gets rank j. The remaining \Vt^ | — (j — 1) vertices in 7^ 
and \Vt^, \ — 1 vertices in Ty have to be shuffled together. According to Remark 
this can be done in 



'\VtJ- [j -1) + \VtJ-1\ flVrJ + lVrJ-j 
different ways. Thus overall there are: 



'[r{u)^^]■\r{%)\■\r{%)\ 



V \VrJ-l 
nVrJ + lVrJ-J 



V |VtJ-1 

different rank functions on T with r(u) = i and r(v) = j. For the probability 

P[r(u) = i, r(v) — j]: 



'[r{u) = i, r(v) = j] 



nu)^^-H%)\-\r{%)\-(y^!j^i^-') 



E,,P[r(^) = ^] ■ \ri%)\ ■ \ri%)\ ■ C'^fJJll^-') 

Since \r{Tu)\ and |r(7^)| are independent of i and j, those factors cancel out, 
giving 

P[r(u)=^].(l^-|+lf-|-^) 
F[r{u) = ^, r{v) = j] = 1^;^-] , , (3) 



Furthermore, note that 

\VtJ + \VtJ-J 
|VrJ-l 



(|Vr|-j)! 



i\VTA-m\VT\-j-i\VrJ-l)y. 

Again, since (|Vr„| — 1)! is independent of i and j, this factor cancels out, and 
so 

pr , , . , , , nr{u) = ^]■U^to~\\yr\-J-k) 

F[r{u) = t, r{v) = j] = '•/ i 9 » 

E,,P[K^)-*]-nLl"o' (l^r|-J-fc) 

Let n = {{ij) -.i < e {1, . . . , \V\}, \VtJ > With this notation, the 

expected edge length E[X] is 

E[X] = E[X\r{u) = i, r{v) = j]F[r{u) = i, r(t;) = j] 

nriu)=^]■U^:o^~'{\VT\-J-k) 



E 

(j,i)Gf2 



T.1 



1 



i + k 



M") = *]-nL="o'"'(i^Ti-j-fc) 





■nriu) = ^]■U!kIt\\Vr\-3-k)_ 




PK^)=*]-n&''"'(|Vr|-J-fc) 





(4) 



9 



Remark 6. Equation 0| enables the estimation of the length of every interior 
edge. For pendant edges, the approach above gives no definite answer. All we 
know is that the time from the latest interior vertex, which has rank n — 1, until 
today is expected to be at most 1/n where n is the number of leaves. 

Suppose that the growth process is stopped as soon as the n — 1-st speciation 
event occurs. In this case the expected length X of a pendant edge below an 
interior vertex v is: 



n— 1 n — 2 

E[X] = Y,F[r{v)^^]Y,— 



k + 1 

l—i. K — l 

The expected depth of vertex v from the first branchpoint is: 

i=l k=l 

So the depth Y of the leaf in question from the first branchpoint has expectation 
independent of v. 

n-2 



n-1 

E 

i=l 


F[r{v) 


i-1 

-IE 

k=l 


n-1 

E: 

1=1 


F[r{v) 


n-2 

-IE 

k=l 


n-2 

E 

fc=i 


1 

k + 1 





k—i 



In other words, assigning to each edge of a given tree topology its expected 
length gives a tree which obeys a molecular clock. 

Remark 7. Often, an inferred tree has vertices with more than two descen- 
dants, i.e. there is lack of resolution due to, e.g. confliciting data. Our calcula- 
tion for the expected edge length assumes a binary tree though. 

However, the expected edge length may be calculated for each possible binary 
resolution of the supertree. Assume the supertree T has the possible binary 
resolutions 71, ... , Tm. For an edge (u, v) in T where u is the immediate ancestor 
of V, the expected edge length is calculated in the trees 7i for i = 1, . . . , to. The 
expected edge length in % is denoted by a for i = 1, . . . ,to. Note that if u 
is a vertex with more than two descendants in T then v is in general not a 
direct descendant of m in 7^. The value in resolution % is then the sum of all 
expected edge lengths on the path from it to w in 7^. 

Calculate the expected edge length E,[X] of {u,v) in the supertree T by 

E[X] - ^^f^^ (5) 
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where the probabiUty of a tree T under the Yule model is [BrownL 1 1994) 

2n-l 



Again, once the expected length of pendant edges is included the resulting 
tree obeys a molecular clock, meaning that all leaves are at the same depth. 



3.2 The coalescent process 

The edge length estimation in the previous section works for the continuous- 
time Yule model. By changing the method above slightly, we get an edge length 
estimation for the coalescent process. In the coalescent setting, we have 

E\X\r(u) = i,r(v) = 7I = V 7 ttt ; r- 

^ ^ ^ ' ' ^ -^J Z-.^ {i + k)ii + k - 1) 

Therefore, the expected edge length for an interior edge {u, v) can be calculated 
by the following modification of Equation 21 



E[X] = 







fe)(i+fe-i)) 


nr{u) = ^]■U^tt'i\yr\-J-k)_ 




P[r{u) = i] 


•nil"o'''(i^vi-j--fc) 





The calculations in Section lTTI and provide exact values for the expected 
length of an interior edge under the Yule or coalescent process as an alterna- 
tive to simulations. However simulations also provide some indication of the 
variability in the estimate of edge lengths, and it may be of interest to also in- 
vestigate analytically the variance (or even the distribution) of the edge length 
in future work, rather than just its mean. 



4 Comparing two interior vertices 

The algorithm RankProb can also be used for comparing two interior vertices. 
Assume again that every rank function on a rooted binary phylogenetic tree T 
is equally likely. The aim is to compare two interior vertices u and v of T. Was 
u more likely before (of lower rank than) v or v before u? In other words, what 
is the probability 

Pn<. nr{u) < r{v)] 

where r{T) is the set of all possible rank functions on T. Note that it does not 
hold P[r(u) < r{v)] — P[r{u) > r{v)] even with the uniform distibution on the 
rank functions. The probability ¥u<v is equivalent to counting all the possible 
rank functions on T in which u has lower rank than v and divide that number 



11 



by all possible rank functions on T. One idea is to sum up the probabilities 
P[r(M) = i,r{v) = j] in Equation O for all i < j which yields to a runtime of 
Odyj"*). The following algorithm Compare solves the problem in quadratic 
time. In the following, for a vertex v, the subtree 7^ of T consists again of v 
and all its descendants. 



Algorithm Compare (T, u, v) 

Input: A rooted binary phylogenetic tree T and two distinct interior vertices 
u and V. 

Output: The probability P„<i, := P[r{u) < r{v)\T]. 
1: Denote the most recent common ancestor of u and v hy pi. 
2: ii pi ^ V then 
3: RETURN P„<„ = 0. 
4: end if 
5: ii pi — u then 
6: RETURN P„<^, = 1. 
7: end if 

8: Let Tp-^ be the subtree of T which is induced by pi. 

9: Delete the vertex pi from 7^^ . The two evolving subtrees are labeled 7^ and 

Ty with u G 7u and v G 
10: Run RankProb(7^, m) and RankProb(7^, w) to get P[r{u) = z] on 7^ and 
¥[r{v) = i] on 7^ for all possible i. 
for i — 1, . . . , \Vt^ I do 

ucum{i) :— Yl\=i '^[t{u) = i] 



11 

12 
13 
14 
15 
16 

17 
18 
19 
20 
21 

22 
23 



end for 



for i ~ 1, . . . , IW,, I do 
for j = 1, . . . \ Vr^ I do 

¥u<v ■= ]Pu<t) + P 

end for 
end for 

tot:= (' "I) 

¥u<v '■— Pu<v/tot 

RETURN P„<„ 



Theorem 8. The algorithm COMPARE returns the value 

F„<„ =P[r(u) <r{v)]. 
The runtime o/ COMPARE is 0{\V\'^). 

Proof. Note that the probability of u having smaller rank than v in tree 7^^ 
equals the probability of u having smaller rank than v in tree T, since for any 
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rank function on 7^^ , there is the same number of hnear extensions to get a 
rank function on the tree T. 

So it is sufficient to calculate the probability Pu<t, in 7^^ . If pi = m then u 
is an ancestor of v in T, so return P„<i, = 1. If pi = t; then v is an ancestor of 
u in T, so return Pu<v = 0. 

Now assume that p\ ^ u and p\ ^ v. The run of RankProb calculates 
the probability P[r(u) = j] in the tree 7^ and P[r(i;) = z] in 7^ for all i. Next, 
combine those two linear orders. Assume that r(v) — i and that j vertices of 
Tu are inserted before v. Inserting j vertices of 7^ into the linear order of 7jj 
before v is possible in (*^^^"') ways (see Remark^. Putting the remaining 

vertices in a linear order is possible in ('^^" 'y^^j'^p' ' ways. The probability 
that the vertex u is among the j vertices which have smaller rank than v is 
P[r(u) < j] = ucura(j). There are |7'(7^)| possible linear orders on 7^ and 
|r(7^)| possible linear orders on Ty. The number of linear orders where vertex 
V has rank i in 7^, u has rank i + j in 7^^ and r(u) < i + j therefore equals 

p^,,. = P[r{v) = ^] ■ \r{%)\ ' (' " J ^ ^) ' ('^''" ' "^^"^i ^^J" ' " ^) -ucumij) ■ \r{T^)\ 

Adding up the p' for each i and j gives the number of linear orders where u has 
smaller rank than v. 

Combining a linear order on 7^ with a linear order on 7^ is possible in 



tot := 



'\VtA + \VtA 
iVrJ 

different ways (see Remark^. There are |r(7^i)| linear orders on 7^ and |?'(7jj)| 
linear orders on 7^ , so on 7^^ , there are 

tot' ('M+|^^"')|r(7;)||r(r.)| 
linear orders. Therefore: 

tot' ~ tot 

with = P[r{v) = i] ■ {'-]+^) ■ {^^^""^^^'^Z^"^'^) ■ ucumij). This shows that 
Compare works correct. 

Since RankProb has quadratic runtime, Compare also has quadratic run- 
time. □ 
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